function [mm_D4,param_D4,policy_D4,perc_college_O,avg_h_O,alpha_O,perc_college_N,avg_h_N,alpha_N,avg_diff_N,avg_diff_O,avg_diff_L,avg_diff_H] = DO_Appendix_D4
    
    %% ==================== Set Parameters and Grids ========================%
    vfi_toler = 0.0001;
    dist_toler = 0.000001;
    
    nu = 0.00622206;              % Quarterly probability of exiting the model
    beta_discount = 0.99;         % Discount factor
    beta = beta_discount*(1-nu);  % Death-adjusted discount factor (factors in exit probability)
    
    alpha(1) = 0.601501;           % Production function parameter - dependent on sector
    alpha(2) = 0.623757;           % Production function parameter - dependent on sector
    d = 0.051475;             % Human capital depreciation rate
    ell = 1.27;                % Matching function parameter
    lambda = 0.246414;        % Probability of search OTJ
    chi = 3.952712;                  % Investment cost is (chi*i)^gamma
    gamma =2;                 % Investment cost is (chi*i)^gamma
    b=0.663351;                   % Unemployment benefit
    muA = 0.104154;               % Staffing agency fee
    
    ns = 2;                   % Number of separation types
    delta = 0.084;         % Separation rate 
    
    c(1) = 0.649976;               % Cost of posting a vacancy - dependent on separation type
    c(2) = 0.744724;
    
    h_mu = 1;                 % Mean of initial human capital distribution
    h_sigma =0.230979;         % Standard deviation of initial human capital distribution
    
    nphi = 2;                 % Number of possible phi values (fixed productivity types)
    pr_phi(1) = 1.0-0.3084251;  % Probability of entering the model with phi=1 (low fixed productivity type)
    pr_phi(2) = 0.3084251;      % Probability of entering the model with phi=2 (high fixed productivity type)
    
    z(1) = 1;                 % Productivity parameter f(phi,h) = z(phi)h^(alpha)
    z(2) = 1.512567;
    
    hlb = 0;               % Human capital grid lower bound 
    hub = 7;               % Human capital grid upper bound 
    nh = 3000;                 % Number of possible human capital values
    h_grid = linspace(hlb,hub,nh); % Create human capital grid
    
    % Need to follow "respondents" from entrance because SS distribution does not save the eta values %
    simul_size = 20000;     % Number of individuals to follow in simulated sample
    nt = 400;               % Number of quarters to follow each individual - 100 yrs max - most will exit the model before this
    simul(simul_size,nt) = struct;
    
    npolicy = 2;
    policy_urate = zeros(npolicy,1);
    policy_GDP = zeros(npolicy,1);
    policy_avg_hc = zeros(npolicy,1);
    policy_avg_wage = zeros(npolicy,1);
    policy_util = zeros(npolicy,1);
    policy_entrant = zeros(npolicy,1);
    
    for policy_ind = 1:npolicy
    
        if (policy_ind==2)
            muA=1;
        end
    
        % Other arrays of use
        xN = zeros(1,ns);
        xO = zeros(1,ns);
        etaO = zeros(1,ns);
        etaN = zeros(1,ns);
        t_Uhat_O = zeros(1,ns);
        t_Uhat_N = zeros(1,ns);
        thetaN = zeros(1,ns);
        pN = zeros(1,ns);
        qN = zeros(1,ns);
        
        xN_N = zeros(1,ns);
        etaN_N = zeros(1,ns);
        t_VNhat_O = zeros(1,ns);
        t_VOhat_O = zeros(1,ns);
        thetaN_N = zeros(1,ns);
        pN_N = zeros(1,ns);
        qN_N = zeros(1,ns);
        t_VNhat_N = zeros(1,ns);
        
        thetaO_N = zeros(1,ns);
        pO_N = zeros(1,ns);
        qO_N = zeros(1,ns);
        t_VOhat_N = zeros(1,ns);
        xO_N = zeros(1,ns);
        etaO_N = zeros(1,ns);
        % ======================= Create discretized h distribution ==================== %
        cdf_h = zeros(1,nh);      % CDF of initial h distribution
        pdf_h = zeros(1,nh);      % PDF of initial h distribution
        
        h_bin_size = (hub-hlb)/(nh-1);
        for ind_h = 1:nh
            if (ind_h ==1)
                bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
                cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
                pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            else
                bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
                bin_ub_last = h_grid(ind_h-1) + (h_bin_size)/2;
                cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
                pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma) - normcdf(bin_ub_last,h_mu,h_sigma);
            end
        end
        %figure
        %plot(h_grid,pdf_h)
        
        %% ============== Backward Loop over Ages to Get Value and Policy Functions ================ %%
        %---- Prespecify Value Functions
        
        U = zeros(nphi,nh);                  % Value of unemployment in production phase                     
        U_hat  = zeros(nphi,nh);             % Value of unemployment in search phase - when searching for optimal job type
        U_hat_update  = zeros(nphi,nh);
        
        VN = zeros(ns,nphi,nh);              % Value of employment in non-outsourced (N) job in production phase
        VO = zeros(ns,nphi,nh);              % Value of employment in outsourced (O) job in production phase
        VN_hat = zeros(ns,nphi,nh);          % Value of employment in N job in search phase
        VO_hat = zeros(ns,nphi,nh);          % Value of employment in O job in search phase
        VN_hat_update = zeros(ns,nphi,nh);
        VO_hat_update = zeros(ns,nphi,nh);
        A_staffing = zeros(ns,nphi,nh);
        %---- Prespecify Value Functions
        opt_ind_hp_N = ones(ns,nphi,nh);     % Optimal index for next period h value (h') if in non-outsourced job
        opt_ind_hp_O = ones(ns,nphi,nh);     % Opitmal index for next period h value (h') if in outsourced job
        
        GU = zeros(nphi,nh);                 % Indicates optimal job type to search for if unemployed - 1=non-outsourced (N) and 2=outsourced (O)
        GN = zeros(ns,nphi,nh);              % Indicates optimal job type to search for if employed in N job - 1=non-outsourced and 2=outsourced
        GO = zeros(ns,nphi,nh);              % Indicates optimal job type to search for if employed in O job - 1=non-outsourced and 2=outsourced
        
        opt_p_U = zeros(nphi,nh);            % Optimal job-finding rate associated with optimal search choice when unemployed
        opt_p_N = zeros(ns,nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in N job
        opt_p_O = zeros(ns,nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in O job
        
        opt_x_U = zeros(nphi,nh);            % x contract value associated with optimal search choice when unemployed
        opt_x_N = zeros(ns,nphi,nh);         % x contract value associated with optimal search choice when employed in N job
        opt_x_O = zeros(ns,nphi,nh);         % x contract value associated with optimal search choice when employed in O job
        
        opt_eta_U = zeros(nphi,nh);          % eta (surplus rate) associated with optimal search choice when unemployed
        opt_eta_N = zeros(ns,nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in N job 
        opt_eta_O = zeros(ns,nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in O job
        
        opt_theta_U = zeros(nphi,nh);        % theta (market tightness) associated with optimal search choice when unemployed
        opt_theta_N = zeros(ns,nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in N job 
        opt_theta_O = zeros(ns,nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in O job
        
        % Find index for next period human capital (h') if the agent only lets human capital depreciate
        ind_hp_min = ones(1,nh);          % Minimum index for h' (index for h' if human capital only depreciates)
        for ind_h = 1:nh
            h = h_grid(ind_h);
            error_min = 100;
            for ind_hp = 1:nh
                hp = h_grid(ind_hp);       % hp is the value of h'
                error = abs(hp - h*(1-d)); % Find index that minimizes (hp-h*(1-d))
                if (error<error_min)
                    error_min = error;
                    ind_hp_min(ind_h) = ind_hp;
                end
            end
        end
        
        %--------------- Begin Value Function Iteration ----------------%
        vfi_error = 1; % Set vfi_error > vfi_toler
        while (vfi_error > vfi_toler)
            for ind_phi = 1:nphi          % Loop over phi types
                for ind_h = 1:nh          % Loop over current h values
                    h = h_grid(ind_h);    % h = current human capital (h) value 
                    
                    b_effect = b*z(ind_phi)*(h^alpha(1));  % Unemployment value
                    c_effect(1) = c(1);                       % Search cost - dependent on separation type
                    c_effect(2) = c(2);
        
                    % ---------- Solve hc Investment Decision ----------%
                    for ind_s = 1:ns % Start loop over possible separation s types 
                        max_VN = -9999;
                        max_VO = -9999;
        
                        for ind_hp = ind_hp_min(ind_h):nh      % Loop over possible next period hp (h') values
                            hp = h_grid(ind_hp);
        
                            temp_VN = z(ind_phi)*(h^alpha(ind_s)) - ((hp - h*(1-d))*chi)^gamma ...
                                + beta*(1-delta)*VN_hat(ind_s,ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp);
        
                            temp_VO = (1-muA)*(z(ind_phi)*(h^alpha(ind_s))) - ((hp - h*(1-d))*chi)^gamma ...
                                + beta*(1-delta)*VO_hat(ind_s,ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp);
        
                            if (temp_VN>max_VN)                % Save ind_hp where N job employment value is highest
                                max_VN = temp_VN;
                                VN(ind_s,ind_phi,ind_h) = temp_VN;
                                opt_ind_hp_N(ind_s,ind_phi,ind_h) = ind_hp;
                            end
                            if (temp_VO>max_VO)                % Save ind_hp where O job employment value is highest
                                max_VO = temp_VO;
                                VO(ind_s,ind_phi,ind_h) = temp_VO;
                                opt_ind_hp_O(ind_s,ind_phi,ind_h) = ind_hp;
                            end
        
                        end % end loop over hp options
                        ind_hp = opt_ind_hp_O(ind_s,ind_phi,ind_h);
                        A_staffing(ind_s,ind_phi,ind_h) = muA*z(ind_phi)*(h^alpha(ind_s)) + (1-delta)*beta*(1-lambda*opt_p_O(ind_s,ind_phi,ind_h))*A_staffing(ind_s,ind_phi,ind_hp);
        
                    end % end loop over possible separation s types
                    U(ind_phi,ind_h) = b_effect + beta*U_hat(ind_phi,ind_hp_min(ind_h));
        
                    %---------------------------------------------------%
                    % ------- Solve Search Decision of Unemployed ----- %
                    for ind_sp = 1:ns % loop over possible future s options
                        xO(ind_sp) = VO(ind_sp,ind_phi,ind_h) - c_effect(ind_sp);
                        etaO(ind_sp) = (xO(ind_sp) - U(ind_phi,ind_h))/(VO(ind_sp,ind_phi,ind_h) - U(ind_phi,ind_h));
                        t_Uhat_O(ind_sp) = xO(ind_sp);
        
                        thetaN(ind_sp) = (((VN(ind_sp,ind_phi,ind_h) - U(ind_phi,ind_h))/c_effect(ind_sp))^(ell/(1+ell)) - 1)^(1/ell);
                        if (~isreal(thetaN(ind_sp)))
                            thetaN(ind_sp)=0;
                            qN(ind_sp) = 0;
                            pN(ind_sp) = 0;
                            xN(ind_sp) = 0;
                            etaN(ind_sp) = 0;
                        else
                            pN(ind_sp) = thetaN(ind_sp)/((1+thetaN(ind_sp)^ell)^(1/ell));
                            qN(ind_sp) = pN(ind_sp)/thetaN(ind_sp);
                            xN(ind_sp) = VN(ind_sp,ind_phi,ind_h) - (c_effect(ind_sp)/qN(ind_sp));
                            etaN(ind_sp) = (xN(ind_sp) - U(ind_phi,ind_h))/(VN(ind_sp,ind_phi,ind_h) - U(ind_phi,ind_h));
                        end
                        t_Uhat_N(ind_sp) = pN(ind_sp)*xN(ind_sp) + (1-pN(ind_sp))*U(ind_phi,ind_h);
        
                    end % end loop over possible future s options
                    
                    % Now Consider which job type's value is highest
                    if (max(t_Uhat_N(:))>max(t_Uhat_O(:))) % If top N job is better than top O job
                        for ind_sp = 1:ns
                            if (t_Uhat_N(ind_sp)==max(t_Uhat_N(:)))
                                GU(ind_phi,ind_h) = 10 + ind_sp;
                                opt_p_U(ind_phi,ind_h) = pN(ind_sp);
                                U_hat_update(ind_phi,ind_h) = t_Uhat_N(ind_sp);
                                opt_x_U(ind_phi,ind_h) = xN(ind_sp);
                                opt_eta_U(ind_phi,ind_h) = etaN(ind_sp);
                                opt_theta_U(ind_phi,ind_h) = thetaN(ind_sp);
                            end
                        end
                    elseif (max(t_Uhat_O(:))>max(t_Uhat_N(:))) % If top O job is better than top N job
                        for ind_sp = 1:ns
                            if (t_Uhat_O(ind_sp)==max(t_Uhat_O(:)))
                                GU(ind_phi,ind_h) = 20 + ind_sp;
                                opt_p_U(ind_phi,ind_h) = 1.0;
                                U_hat_update(ind_phi,ind_h) = t_Uhat_O(ind_sp);
                                opt_x_U(ind_phi,ind_h) = xO(ind_sp);
                                opt_eta_U(ind_phi,ind_h) = etaO(ind_sp);
                                opt_theta_U(ind_phi,ind_h) = 1.0;
                            end
                        end
                    else % If top O and N jobs are the same
                        for ind_sp = 1:ns
                            if (t_Uhat_N(ind_sp)==max(t_Uhat_N(:)))
                                GU(ind_phi,ind_h) = 30 + ind_sp;
                                opt_p_U(ind_phi,ind_h) = pN(ind_sp);
                                U_hat_update(ind_phi,ind_h) = t_Uhat_N(ind_sp);
                                opt_x_U(ind_phi,ind_h) = xN(ind_sp);
                                opt_eta_U(ind_phi,ind_h) = etaN(ind_sp);
                                opt_theta_U(ind_phi,ind_h) = thetaN(ind_sp);
                            end
                        end
                    end
        
                    %---------------------------------------------------%
                    % ------- Solve Search Decision of Employed --------%
                    for ind_s = 1:ns % loop over current s value
                        for ind_sp = 1:ns
        
                            xO(ind_sp) = VO(ind_sp,ind_phi,ind_h) - c_effect(ind_sp);
                            etaO(ind_sp) = (xO(ind_sp) - U(ind_phi,ind_h))/(VO(ind_sp,ind_phi,ind_h) - U(ind_phi,ind_h));
        
                            t_VNhat_O(ind_sp) = lambda*xO(ind_sp) + (1-lambda)*VN(ind_s,ind_phi,ind_h);
                            t_VOhat_O(ind_sp) = lambda*xO(ind_sp) + (1-lambda)*VO(ind_s,ind_phi,ind_h);
        
                            thetaN_N(ind_sp) = (((VN(ind_sp,ind_phi,ind_h) - VN(ind_s,ind_phi,ind_h))/c_effect(ind_sp))^(ell/(1+ell)) - 1)^(1/ell);
                            if (~isreal(thetaN_N(ind_sp)))
                                thetaN_N(ind_sp) =0;
                                qN_N(ind_sp) = 0;
                                pN_N(ind_sp) = 0;
                                xN_N(ind_sp) = 0;
                                etaN_N(ind_sp) = 0;
                            else
                                pN_N(ind_sp) = thetaN_N(ind_sp)/((1+thetaN_N(ind_sp)^ell)^(1/ell));
                                qN_N(ind_sp) = pN_N(ind_sp)/thetaN_N(ind_sp);
                                xN_N(ind_sp) = VN(ind_sp,ind_phi,ind_h) - (c_effect(ind_sp)/qN_N(ind_sp));
                                etaN_N(ind_sp) = (xN_N(ind_sp) - U(ind_phi,ind_h))/(VN(ind_sp,ind_phi,ind_h) - U(ind_phi,ind_h));
                            end
                            t_VNhat_N(ind_sp) = lambda*pN_N(ind_sp)*xN_N(ind_sp) + (1-lambda*pN_N(ind_sp))*VN(ind_s,ind_phi,ind_h);
        
                            thetaO_N(ind_sp) = (((VN(ind_sp,ind_phi,ind_h) - VO(ind_s,ind_phi,ind_h))/c_effect(ind_sp))^(ell/(1+ell)) - 1)^(1/ell);
                            if (~isreal(thetaO_N(ind_sp)))
                                thetaO_N(ind_sp)=0;
                                qO_N(ind_sp) = 0;
                                pO_N(ind_sp) = 0;
                                xO_N(ind_sp) = 0;
                                etaO_N(ind_sp) = 0;
                            else
                                pO_N(ind_sp) = thetaO_N(ind_sp)/((1+thetaO_N(ind_sp)^ell)^(1/ell));
                                qO_N(ind_sp) = pO_N(ind_sp)/thetaO_N(ind_sp);
                                xO_N(ind_sp) = VN(ind_sp,ind_phi,ind_h) - (c_effect(ind_sp)/qO_N(ind_sp));
                                etaO_N(ind_sp) = (xO_N(ind_sp) - U(ind_phi,ind_h))/(VN(ind_sp,ind_phi,ind_h) - U(ind_phi,ind_h));
                            end  
                            t_VOhat_N(ind_sp) = lambda*pO_N(ind_sp)*xO_N(ind_sp) + (1-lambda*pO_N(ind_sp))*VO(ind_s,ind_phi,ind_h);
        
                        end % end loop over possible future s values
                
        
                        % For those Currently in N - Now Consider which job type is highest
                        if (max(t_VNhat_N(:))>max(t_VNhat_O(:))) % If top N job is better than top O job
                            for ind_sp = 1:ns
                                if (t_VNhat_N(ind_sp)==max(t_VNhat_N(:)))
                                    GN(ind_s,ind_phi,ind_h) = 10 + ind_sp;
                                    opt_p_N(ind_s,ind_phi,ind_h) = pN_N(ind_sp);
                                    VN_hat_update(ind_s,ind_phi,ind_h) = t_VNhat_N(ind_sp);
                                    opt_x_N(ind_s,ind_phi,ind_h) = xN_N(ind_sp);
                                    opt_eta_N(ind_s,ind_phi,ind_h) = etaN_N(ind_sp);
                                    opt_theta_N(ind_s,ind_phi,ind_h) = thetaN_N(ind_sp);
                                end
                            end
                        elseif (max(t_VNhat_O(:))>max(t_VNhat_N(:))) % If top O job is better than top N job
                            for ind_sp = 1:ns
                                if (t_VNhat_O(ind_sp)==max(t_VNhat_O(:)))
                                    GN(ind_s,ind_phi,ind_h) = 20 + ind_sp;
                                    opt_p_N(ind_s,ind_phi,ind_h) = 1.0;
                                    VN_hat_update(ind_s,ind_phi,ind_h) = t_VNhat_O(ind_sp);
                                    opt_x_N(ind_s,ind_phi,ind_h) = xO(ind_sp);
                                    opt_eta_N(ind_s,ind_phi,ind_h) = etaO(ind_sp);
                                    opt_theta_N(ind_s,ind_phi,ind_h) = 1.0;
                                end
                            end
                        else % If top O and N jobs are the same
                            for ind_sp = 1:ns
                                if (t_VNhat_N(ind_sp)==max(t_VNhat_N(:)))
                                    GN(ind_s,ind_phi,ind_h) = 30 + ind_sp;
                                    opt_p_N(ind_s,ind_phi,ind_h) = pN_N(ind_sp);
                                    VN_hat_update(ind_s,ind_phi,ind_h) = t_VNhat_N(ind_sp);
                                    opt_x_N(ind_s,ind_phi,ind_h) = xN_N(ind_sp);
                                    opt_eta_N(ind_s,ind_phi,ind_h) = etaN_N(ind_sp);
                                    opt_theta_N(ind_s,ind_phi,ind_h) = thetaN_N(ind_sp);
                                end
                            end
                        end
        
                        % For those Currently in O - Now Consider which job type is highest
                        if (max(t_VOhat_N(:))>max(t_VOhat_O(:))) % If top N job is better than top O job
                            for ind_sp = 1:ns
                                if (t_VOhat_N(ind_sp)==max(t_VOhat_N(:)))
                                    GO(ind_s,ind_phi,ind_h) = 10 + ind_sp;
                                    opt_p_O(ind_s,ind_phi,ind_h) = pO_N(ind_sp);
                                    VO_hat_update(ind_s,ind_phi,ind_h) = t_VOhat_N(ind_sp);
                                    opt_x_O(ind_s,ind_phi,ind_h) = xO_N(ind_sp);
                                    opt_eta_O(ind_s,ind_phi,ind_h) = etaO_N(ind_sp);
                                    opt_theta_O(ind_s,ind_phi,ind_h) = thetaO_N(ind_sp);
                                end
                            end
                        elseif (max(t_VOhat_O(:))>max(t_VOhat_N(:))) % If top O job is better than top N job
                            for ind_sp = 1:ns
                                if (t_VOhat_O(ind_sp)==max(t_VOhat_O(:)))
                                    GO(ind_s,ind_phi,ind_h) = 20 + ind_sp;
                                    opt_p_O(ind_s,ind_phi,ind_h) = 1.0;
                                    VO_hat_update(ind_s,ind_phi,ind_h) = t_VOhat_O(ind_sp);
                                    opt_x_O(ind_s,ind_phi,ind_h) = xO(ind_sp);
                                    opt_eta_O(ind_s,ind_phi,ind_h) = etaO(ind_sp);
                                    opt_theta_O(ind_s,ind_phi,ind_h) = 1.0;
                                end
                            end
                        else % If top O and N jobs are the same
                            for ind_sp = 1:ns
                                if (t_VOhat_N(ind_sp)==max(t_VOhat_N(:)))
                                    GO(ind_s,ind_phi,ind_h) = 30 + ind_sp;
                                    opt_p_O(ind_s,ind_phi,ind_h) = pO_N(ind_sp);
                                    VO_hat_update(ind_s,ind_phi,ind_h) = t_VOhat_N(ind_sp);
                                    opt_x_O(ind_s,ind_phi,ind_h) = xO_N(ind_sp);
                                    opt_eta_O(ind_s,ind_phi,ind_h) = etaO_N(ind_sp);
                                    opt_theta_O(ind_s,ind_phi,ind_h) = thetaO_N(ind_sp);
                                end
                            end
                        end
        
                    end % end loop over current s value
        
                end % End loop over h values
            end % End loop over phi values
        
            error_U_hat = max(max(abs(U_hat_update(:,:) - U_hat(:,:))))/max(max(abs(U_hat_update(:,:))));
            error_VN_hat = max(max(max(abs(VN_hat_update(:,:,:) - VN_hat(:,:,:)))))/max(max(max(abs(VN_hat_update(:,:,:)))));
            error_VO_hat = max(max(max(abs(VO_hat_update(:,:,:) - VO_hat(:,:,:)))))/max(max(max(abs(VO_hat_update(:,:,:)))));
        
            vfi_error = max([error_U_hat,error_VN_hat,error_VO_hat]);
        
            U_hat(:,:) = U_hat_update(:,:);
            VN_hat(:,:,:) = VN_hat_update(:,:,:);
            VO_hat(:,:,:) = VO_hat_update(:,:,:);
        
        end % End while (vfi_error > vfi_toler)  
        
        %==================================================================================================================%
        %% ================================ Solve for Steady State Distribution ========================================= %%
        %==================================================================================================================%
            
        dist_U = zeros(nphi,nh);              % Steady state mass of unemployed agents
        dist_N = zeros(ns,nphi,nh);           % Steady state mass of agents employed in non-outsourced job
        dist_O = zeros(ns,nphi,nh);           % Steady state mass of agents employed in outsourced job
        
        dist_U_update = zeros(nphi,nh);
        dist_N_update = zeros(ns,nphi,nh);
        dist_O_update = zeros(ns,nphi,nh);
        
        dist_U(1,1) = pr_phi(1);           % Start with a guess regarding how the mass of agents is distributed
        dist_U(2,1) = pr_phi(2);
        
        dist_error = 1;
        while (dist_error>dist_toler)      % Iterate until distribution guess equals resulting distribution
        
            for ind_phi = 1:nphi           % Loop over phi types
                for ind_h = 1:nh           % Loop over current h values 
                    %=================================%
                    %========= New Entrants ==========%
                    %=================================%
                    % New Entrant - Finds Job
                    if (GU(ind_phi,ind_h)<20 || GU(ind_phi,ind_h)>=30)     % Searches for N
                        for ind_sp = 1:ns
                            if (ind_sp==(GU(ind_phi,ind_h)-10) || ind_sp==(GU(ind_phi,ind_h)-30))
                                dist_N_update(ind_sp,ind_phi,ind_h) = dist_N_update(ind_sp,ind_phi,ind_h)...
                                    + opt_p_U(ind_phi,ind_h)*pr_phi(ind_phi)*nu*pdf_h(ind_h);
                            end
                        end
                    else                                                   % Searches for O
                        for ind_sp = 1:ns
                            if (ind_sp==(GU(ind_phi,ind_h)-20) )
                                dist_O_update(ind_sp,ind_phi,ind_h) = dist_O_update(ind_sp,ind_phi,ind_h)...
                                    + opt_p_U(ind_phi,ind_h)*pr_phi(ind_phi)*nu*pdf_h(ind_h);
                            end
                        end
                    end
                    % New Entrant - Does Not Find Job
                    dist_U_update(ind_phi,ind_h) = dist_U_update(ind_phi,ind_h) + (1-opt_p_U(ind_phi,ind_h))*pr_phi(ind_phi)*nu*pdf_h(ind_h);
        
                    %=================================%
                    %===== Previously Unemployed =====%
                    %=================================%
                    ind_hp = ind_hp_min(ind_h);             % New h value for the unemployed is what they get from allowing human capital to depreciate
        
                    % Unemployed - Finds a job
                    if (GU(ind_phi,ind_hp)<20 || GU(ind_phi,ind_hp)>=30)     % Searches for N
                        for ind_sp = 1:ns
                            if (ind_sp==(GU(ind_phi,ind_hp)-10) || ind_sp==(GU(ind_phi,ind_hp)-30))
                                dist_N_update(ind_sp,ind_phi,ind_hp) = dist_N_update(ind_sp,ind_phi,ind_hp)...
                                    + opt_p_U(ind_phi,ind_hp)*dist_U(ind_phi,ind_h)*(1-nu);
                            end
                        end
                    else                    % Searches for O
                        for ind_sp = 1:ns 
                            if (ind_sp==(GU(ind_phi,ind_hp)-20) )
                                dist_O_update(ind_sp,ind_phi,ind_hp) = dist_O_update(ind_sp,ind_phi,ind_hp)...
                                    + opt_p_U(ind_phi,ind_hp)*dist_U(ind_phi,ind_h)*(1-nu);
                            end
                        end
                    end 
                    % Unemployed - Does not find a job
                    dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp) + (1-opt_p_U(ind_phi,ind_hp))*dist_U(ind_phi,ind_h)*(1-nu);
        
                    for ind_s = 1:ns % For employed: loop over s in their last period job
                        %=================================%
                        %======= Previously in N =========%
                        %=================================%
                        ind_hp = opt_ind_hp_N(ind_s,ind_phi,ind_h);   % New h value for employed in N job is the result of their optimal investment choice
        
                        % ----- N - No Separation Shock - Finds New Job
                        if (GN(ind_s,ind_phi,ind_hp)<20 || GN(ind_s,ind_phi,ind_hp)>=30)      % N - Search OTJ for N
                            for ind_sp = 1:ns
                                if (ind_sp==(GN(ind_s,ind_phi,ind_hp)-10) || ind_sp==(GN(ind_s,ind_phi,ind_hp)-30))
                                    dist_N_update(ind_sp,ind_phi,ind_hp) = dist_N_update(ind_sp,ind_phi,ind_hp)...
                                        + lambda*opt_p_N(ind_s,ind_phi,ind_hp)*(1-delta)*dist_N(ind_s,ind_phi,ind_h)*(1-nu);
                                end
                            end 
                        else                                                                  % N - Search OTJ for O
                            for ind_sp = 1:ns
                                if (ind_sp==(GN(ind_s,ind_phi,ind_hp)-20))
                                    dist_O_update(ind_sp,ind_phi,ind_hp) = dist_O_update(ind_sp,ind_phi,ind_hp)...
                                        + lambda*opt_p_N(ind_s,ind_phi,ind_hp)*(1-delta)*dist_N(ind_s,ind_phi,ind_h)*(1-nu);
                                end
                            end
                        end
        
                        % ----- N - No Separation Shock - Does Not Find New Job
                        dist_N_update(ind_s,ind_phi,ind_hp) = dist_N_update(ind_s,ind_phi,ind_hp)...
                            + (1-lambda*opt_p_N(ind_s,ind_phi,ind_hp))*(1-delta)*dist_N(ind_s,ind_phi,ind_h)*(1-nu);
        
                        % ----- N - Separation Shock - Finds Job
                        if (GU(ind_phi,ind_hp)<20 || GU(ind_phi,ind_hp)>=30)      % N - Separation - Search for N
                            for ind_sp = 1:ns
                                if (ind_sp==(GU(ind_phi,ind_hp)-10) || ind_sp==(GU(ind_phi,ind_hp)-30))
                                    dist_N_update(ind_sp,ind_phi,ind_hp) = dist_N_update(ind_sp,ind_phi,ind_hp)...
                                        + opt_p_U(ind_phi,ind_hp)*delta*dist_N(ind_s,ind_phi,ind_h)*(1-nu);
                                end 
                            end
                        else                                                       % N - Separation - Search for O
                            for ind_sp = 1:ns
                                if (ind_sp==(GU(ind_phi,ind_hp)-20) )
                                    dist_O_update(ind_sp,ind_phi,ind_hp) = dist_O_update(ind_sp,ind_phi,ind_hp)...
                                        + opt_p_U(ind_phi,ind_hp)*delta*dist_N(ind_s,ind_phi,ind_h)*(1-nu);
                                end
                            end
                        end
        
                        % ----- N - Separation Shock - Does Not Find New Job
                        dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp)...
                            + (1-opt_p_U(ind_phi,ind_hp))*delta*dist_N(ind_s,ind_phi,ind_h)*(1-nu);
        
                        %=================================%
                        %======= Previously in O =========%
                        %=================================%
                        ind_hp = opt_ind_hp_O(ind_s,ind_phi,ind_h);   % New h value for employed in O job is the result of their optimal investment choice
        
                        % ----- O - No Separation Shock - Finds New Job
                        if (GO(ind_s,ind_phi,ind_hp)<20 || GO(ind_s,ind_phi,ind_hp)>=30)      % O - Search OTJ for N
                            for ind_sp = 1:ns
                                if (ind_sp==(GO(ind_s,ind_phi,ind_hp)-10) || ind_sp==(GO(ind_s,ind_phi,ind_hp)-30))
                                    dist_N_update(ind_sp,ind_phi,ind_hp) = dist_N_update(ind_sp,ind_phi,ind_hp)...
                                        + lambda*opt_p_O(ind_s,ind_phi,ind_hp)*(1-delta)*dist_O(ind_s,ind_phi,ind_h)*(1-nu);
                                end
                            end
                        else                                                                  % O - Search OTJ for O
                            for ind_sp = 1:ns
                                if (ind_sp==(GO(ind_s,ind_phi,ind_hp)-20))
                                    dist_O_update(ind_sp,ind_phi,ind_hp) = dist_O_update(ind_sp,ind_phi,ind_hp)...
                                        + lambda*opt_p_O(ind_s,ind_phi,ind_hp)*(1-delta)*dist_O(ind_s,ind_phi,ind_h)*(1-nu);
                                end
                            end
                        end
        
                        % ----- O - No Separation Shock - Does Not Find New Job
                        dist_O_update(ind_s,ind_phi,ind_hp) = dist_O_update(ind_s,ind_phi,ind_hp)...
                            + (1-lambda*opt_p_O(ind_s,ind_phi,ind_hp))*(1-delta)*dist_O(ind_s,ind_phi,ind_h)*(1-nu);
        
                        % ----- O - Separation Shock - Finds Job
                        if (GU(ind_phi,ind_hp)<20 || GU(ind_phi,ind_hp)>=30)       % O - Separation - Search for N
                            for ind_sp = 1:ns
                                if (ind_sp==(GU(ind_phi,ind_hp)-10) || ind_sp==(GU(ind_phi,ind_hp)-30))
                                    dist_N_update(ind_sp,ind_phi,ind_hp) = dist_N_update(ind_sp,ind_phi,ind_hp)...
                                        + opt_p_U(ind_phi,ind_hp)*delta*dist_O(ind_s,ind_phi,ind_h)*(1-nu);
                                end
                            end
                        else                                                       % O - Separation - Search for O
                            for ind_sp = 1:ns
                                if (ind_sp==(GU(ind_phi,ind_hp)-20))
                                    dist_O_update(ind_sp,ind_phi,ind_hp) = dist_O_update(ind_sp,ind_phi,ind_hp)...
                                        + opt_p_U(ind_phi,ind_hp)*delta*dist_O(ind_s,ind_phi,ind_h)*(1-nu);
                                end
                            end
                        end
        
                        % ----- O - Separation Shock - Does Not Find New Job
                        dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp)...
                            + (1-opt_p_U(ind_phi,ind_hp))*delta*dist_O(ind_s,ind_phi,ind_h)*(1-nu);
        
                    end % end loop over s in their last period job
        
                 end % End loop over last period h value
            end % End loop over phi values
        
            error_dist_U = max(max(abs(dist_U_update(:,:) - dist_U(:,:))));
            error_dist_N = max(max(max(abs(dist_N_update(:,:,:) - dist_N(:,:,:)))));
            error_dist_O = max(max(max(abs(dist_O_update(:,:,:) - dist_O(:,:,:)))));
        
            dist_error = max([error_dist_U,error_dist_N,error_dist_O]);
        
            dist_U(:,:) = dist_U_update(:,:);
            dist_N(:,:,:) = dist_N_update(:,:,:);
            dist_O(:,:,:) = dist_O_update(:,:,:);
        
            dist_U_update(:,:) = 0;
            dist_N_update(:,:,:) = 0;
            dist_O_update(:,:,:) = 0;
        
        end % end while (dist_error>dist_toler)           
        
        %% =============================================================================================================== %%
        %======================================== Simulation for Wage Moments and Results ==================================%
        %================================================================================================================= %%
        
        rng(123456789);         % Set random number seed so results are exactly the same each time the code is ran     
    
        for id = 1:simul_size   % Loop over each individual in the sample
            for t = 1:nt        % Loop over each time period the individual is observed
        
                if (t==1)
                    % For new entrants, assign phi and h, then see if they are successful in their job search from unemployment %
                    % Assign starting h value
                    r1 = rand;
                    simul(id,t).h_ind = 1;  % If r1 is not bigger than any CDF values, assign to the lowest slot
                    for ind_h = 2:nh
                        if (r1>cdf_h(ind_h-1)) 
                            simul(id,t).h_ind = ind_h;
                        end 
                    end
                    ind_h = simul(id,t).h_ind;
                    h = h_grid(ind_h);
                    % Now assign phi type
                    r1 = rand;
                    if r1 < pr_phi(1)
                        simul(id,t).phi_ind = 1;
                    else
                        simul(id,t).phi_ind = 2;
                    end
                    ind_phi = simul(id,t).phi_ind;
                    % See if new entrant is successful in job search from unemployment - - If successful, save their wage %
                    r1 = rand;
                    if (r1 < opt_p_U(ind_phi,ind_h)) % Finds job
                        if (GU(ind_phi,ind_h)<20 || GU(ind_phi,ind_h)>=30)         % Search for N
                            for ind_sp = 1:ns
                                if (ind_sp == (GU(ind_phi,ind_h)-10) || ind_sp == (GU(ind_phi,ind_h)-30))
                                    simul(id,t).stat = 1;                          % Status "stat" 1 = in non-outsourced job
                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                    eta = simul(id,t).eta;
                                    simul(id,t).s_ind = ind_sp;
                                    ind_s = simul(id,t).s_ind;
        
                                    ind_hp = opt_ind_hp_N(ind_s,ind_phi,ind_h);    % This is the h' value they will choose for next period 
                                    hp = h_grid(ind_hp);
                                    simul(id,t).wage = eta*(VN(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_s,ind_phi,ind_hp)*opt_x_N(ind_s,ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_s,ind_phi,ind_hp))*(eta*VN(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                end
                            end
                        else                                                       % Search for O
                            for ind_sp = 1:ns
                                if (ind_sp == (GU(ind_phi,ind_h)-20))
                                    simul(id,t).stat = 2;                          % Status "stat" 2 = in outsourced job 
                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                    eta = simul(id,t).eta;
                                    simul(id,t).s_ind = ind_sp;
                                    ind_s = simul(id,t).s_ind;
        
                                    ind_hp = opt_ind_hp_O(ind_s,ind_phi,ind_h);    % This is the h' value they will choose for next period 
                                    hp = h_grid(ind_hp);
        
                                    simul(id,t).wage = eta*(VO(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_s,ind_phi,ind_hp)*opt_x_O(ind_s,ind_phi,ind_hp)...
                                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_s,ind_phi,ind_hp))*(eta*VO(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                end
                            end
                        end 
                    else  % Does not find job
                        simul(id,t).stat = 0;                                    % Status "stat" 0 = unemployed  
                        simul(id,t).wage = 0.0;
                        simul(id,t).eta = 0.0;
                        simul(id,t).s_ind=0;
                    end   % End if t=1 individual finds job or not    
        
                else % so t>1
                    simul(id,t).phi_ind = simul(id,t-1).phi_ind;                 % Individual's phi value is permanent
                    ind_phi = simul(id,t).phi_ind;
        
                    if (simul(id,t-1).stat==999) % if dead                       % If the agent already exited the model they remain dead             
                        simul(id,t).stat=999;                                    % Status "stat" 999 = dead
                    else   % was not dead in last period
                        rnu = rand;
                        if (rnu<nu) % Dies at end of last period
                            simul(id,t).stat=999;
                        else        % Does NOT die at end of last period
        
                            if ((simul(id,t-1).stat)==0)        % If unemployed last period 
                                % Recall h from last period and apply any depreciation
                                ind_h = simul(id,t-1).h_ind;
                                simul(id,t).h_ind =  ind_hp_min(ind_h);
                                ind_h = simul(id,t).h_ind;
                                h = h_grid(ind_h);
        
                                % If unemployed last period - see if matches from unemployment !
                                r1 = rand;
                                if (r1 < opt_p_U(ind_phi,ind_h))    % Finds job
                                    if (GU(ind_phi,ind_h)<20 || GU(ind_phi,ind_h)>=30)      % Search for N job
                                        for ind_sp = 1:ns
                                            if (ind_sp == (GU(ind_phi,ind_h) -10) || ind_sp == (GU(ind_phi,ind_h) -30))
                                                simul(id,t).stat = 1; 
                                                simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                                eta = simul(id,t).eta;
                                                simul(id,t).s_ind = ind_sp;
                                                ind_s = simul(id,t).s_ind;
                                                ind_hp = opt_ind_hp_N(ind_s,ind_phi,ind_h);
                                                hp = h_grid(ind_hp);
                                                simul(id,t).wage = eta*(VN(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                    - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_s,ind_phi,ind_hp)*opt_x_N(ind_s,ind_phi,ind_hp)...
                                                    - beta*(1-delta)*(1-lambda*opt_p_N(ind_s,ind_phi,ind_hp))*(eta*VN(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                            end
                                        end
                                    else                                                    % Search for O job
                                        for ind_sp = 1:ns
                                            if (ind_sp == (GU(ind_phi,ind_h) -20))
                                                simul(id,t).stat = 2; 
                                                simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                                eta = simul(id,t).eta;
                                                simul(id,t).s_ind = ind_sp;
                                                ind_s = simul(id,t).s_ind;
        
                                                ind_hp = opt_ind_hp_O(ind_s,ind_phi,ind_h);
                                                hp = h_grid(ind_hp);
                                                simul(id,t).wage = eta*(VO(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                    - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_s,ind_phi,ind_hp)*opt_x_O(ind_s,ind_phi,ind_hp)...
                                                    - beta*(1-delta)*(1-lambda*opt_p_O(ind_s,ind_phi,ind_hp))*(eta*VO(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                            end
                                        end
                                    end
                                else                                 % Unemployed last period and does not find job%
                                    simul(id,t).stat = 0;
                                    simul(id,t).wage = 0.0;
                                    simul(id,t).eta = 0.0;
                                    simul(id,t).s_ind = 0;
                                end
        
                            elseif ((simul(id,t-1).stat)==1)    % If in N last period  
                                ind_s = simul(id,t-1).s_ind;
                                simul(id,t).s_ind = ind_s;
                                % Recall h from last period and apply any investment
                                ind_h = simul(id,t-1).h_ind;
                                simul(id,t).h_ind = opt_ind_hp_N(ind_s,ind_phi,ind_h);
                                ind_h = simul(id,t).h_ind;
                                h = h_grid(ind_h);
        
                                % If N last period - see if separation shock %
                                rdelta = rand;
                                if (rdelta < delta)  % N last period - separation shock    %
                                    % See if finds job from unemployment
                                    r1 = rand;
                                    if (r1 < opt_p_U(ind_phi,ind_h)) % N last period - separation shock - finds job from unemployment
                                        if (GU(ind_phi,ind_h)<20 || GU(ind_phi,ind_h)>=30)      % Search for N
                                            for ind_sp=1:ns
                                                if (ind_sp == (GU(ind_phi,ind_h)-10) || ind_sp == (GU(ind_phi,ind_h)-30))
                                                    simul(id,t).stat = 1; 
                                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                                    eta = simul(id,t).eta;
                                                    simul(id,t).s_ind = ind_sp;
                                                    ind_s = simul(id,t).s_ind;
        
                                                    ind_hp = opt_ind_hp_N(ind_s,ind_phi,ind_h);
                                                    hp = h_grid(ind_hp);
                                                    simul(id,t).wage = eta*(VN(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_s,ind_phi,ind_hp)*opt_x_N(ind_s,ind_phi,ind_hp)...
                                                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_s,ind_phi,ind_hp))*(eta*VN(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                                end
                                            end
                                        else                                                    % Search for O
                                            for ind_sp = 1:ns
                                                if (ind_sp == (GU(ind_phi,ind_h)-20))
                                                    simul(id,t).stat = 2; 
                                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                                    eta = simul(id,t).eta;
        
                                                    simul(id,t).s_ind = ind_sp;
                                                    ind_s = simul(id,t).s_ind;
                                                    ind_hp = opt_ind_hp_O(ind_s,ind_phi,ind_h);
                                                    hp = h_grid(ind_hp);
                                                    simul(id,t).wage = eta*(VO(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_s,ind_phi,ind_hp)*opt_x_O(ind_s,ind_phi,ind_hp)...
                                                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_s,ind_phi,ind_hp))*(eta*VO(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                                end
                                            end
                                        end
                                    else                                 % N last period - separation shock - does find NOT job from unemployment
                                        simul(id,t).stat = 0;
                                        simul(id,t).wage = 0.0;
                                        simul(id,t).eta = 0.0;
                                        simul(id,t).s_ind = 0;
                                    end
                                else                   % N last period - no separation shock %
                                    % See if successfully searchs OTJ
                                    r1 = rand;
                                    if (r1 < lambda*opt_p_N(ind_s,ind_phi,ind_h)) % N last period - no separation shock - successfully searchs OTJ
                                
                                        if (GN(ind_s,ind_phi,ind_h)<20 || GN(ind_s,ind_phi,ind_h)>=30)             % N - Search OTJ for N
                                            for ind_sp = 1:ns
                                                if (ind_sp== (GN(ind_s,ind_phi,ind_h)-10) || ind_sp== (GN(ind_s,ind_phi,ind_h)-30) )
                                                    simul(id,t).stat = 1; 
                                                    simul(id,t).eta = opt_eta_N(ind_s,ind_phi,ind_h);
                                                    eta = simul(id,t).eta;
        
                                                    simul(id,t).s_ind = ind_sp;
                                                    ind_s = simul(id,t).s_ind;
                                                    ind_hp = opt_ind_hp_N(ind_s,ind_phi,ind_h);
                                                    hp = h_grid(ind_hp);
                                                    simul(id,t).wage = eta*(VN(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_s,ind_phi,ind_hp)*opt_x_N(ind_s,ind_phi,ind_hp)...
                                                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_s,ind_phi,ind_hp))*(eta*VN(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                                end
                                            end
                                        else                                                           % N - Search OTJ for O
                                            for ind_sp = 1:ns
                                                if (ind_sp== (GN(ind_s,ind_phi,ind_h)-20))
                                                    simul(id,t).stat = 2; 
                                                    simul(id,t).eta = opt_eta_N(ind_s,ind_phi,ind_h);
                                                    eta = simul(id,t).eta;
        
                                                    simul(id,t).s_ind = ind_sp;
                                                    ind_s = simul(id,t).s_ind;
                                                    ind_hp = opt_ind_hp_O(ind_s,ind_phi,ind_h);
                                                    hp = h_grid(ind_hp);
                                                    simul(id,t).wage = eta*(VO(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_s,ind_phi,ind_hp)*opt_x_O(ind_s,ind_phi,ind_hp)...
                                                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_s,ind_phi,ind_hp))*(eta*VO(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                                end
                                            end
                                        end
                                    else                                         % N last period - no separation shock - does NOT successfully search OTJ 
                                        simul(id,t).stat = simul(id,t-1).stat;
                                        simul(id,t).eta = simul(id,t-1).eta;
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_s,ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_s,ind_phi,ind_hp)*opt_x_N(ind_s,ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_s,ind_phi,ind_hp))*(eta*VN(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                    end
                                end % End if receives separation shock when N last period 
                            else                                % If in O last period 
                                ind_s = simul(id,t-1).s_ind;
                                simul(id,t).s_ind = ind_s;
        
                                % Recall h from last period and apply any investment
                                ind_h = simul(id,t-1).h_ind; 
                                simul(id,t).h_ind = opt_ind_hp_O(ind_s,ind_phi,ind_h);
                                ind_h = simul(id,t).h_ind;
                                h = h_grid(ind_h);
        
                                % If O last period - see if separation shock %
                                rdelta = rand;
                                if (rdelta < delta) % O last period - separation shock 
                                    % See if finds job from unemployment
                                    r1 = rand;
                                    if (r1 < opt_p_U(ind_phi,ind_h))    % N last period - separation shock - finds job from unemployment
                                        if (GU(ind_phi,ind_h)<20 || GU(ind_phi,ind_h)>=30)       % Search for N
                                            for ind_sp = 1:ns
                                                if (ind_sp== (GU(ind_phi,ind_h)-10) || ind_sp== (GU(ind_phi,ind_h)-30))
                                                    simul(id,t).stat = 1; 
                                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                                    eta = simul(id,t).eta;
        
                                                    simul(id,t).s_ind = ind_sp;
                                                    ind_s = simul(id,t).s_ind;
                                                    ind_hp = opt_ind_hp_N(ind_s,ind_phi,ind_h);
                                                    hp = h_grid(ind_hp);
                                                    simul(id,t).wage = eta*(VN(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_s,ind_phi,ind_hp)*opt_x_N(ind_s,ind_phi,ind_hp)...
                                                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_s,ind_phi,ind_hp))*(eta*VN(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                                                end
                                            end
                                        else                                                     % Search for O
                                            for ind_sp = 1:ns
                                                if (ind_sp == (GU(ind_phi,ind_h)-20))
                                                    simul(id,t).stat = 2; 
                                                    simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                                    eta = simul(id,t).eta;
        
                                                    simul(id,t).s_ind = ind_sp;
                                                    ind_s = simul(id,t).s_ind;
                                                    ind_hp = opt_ind_hp_O(ind_s,ind_phi,ind_h);
                                                    hp = h_grid(ind_hp);
                                                    simul(id,t).wage = eta*(VO(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_s,ind_phi,ind_hp)*opt_x_O(ind_s,ind_phi,ind_hp)...
                                                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_s,ind_phi,ind_hp))*(eta*VO(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                                end
                                            end
                                        end
                           
                                    else                                 % O last period - separation shock - does find NOT job from unemployment
                                        simul(id,t).stat = 0;
                                        simul(id,t).wage = 0.0;
                                        simul(id,t).eta = 0.0;
                                        simul(id,t).s_ind = ind_s;
                                    end
                                else                   % O last period - no separation shock %
                                    %See if successfully searchs OTJ
                                    r1 = rand;
                                    if (r1 < lambda*opt_p_O(ind_s,ind_phi,ind_h))   % O last period - no separation shock - successfully searchs OTJ
                                        if (GO(ind_s,ind_phi,ind_h)<20 || GO(ind_s,ind_phi,ind_h)>=30)            % O - Search OTJ for N
                                            for ind_sp = 1:ns
                                                if (ind_sp == (GO(ind_s,ind_phi,ind_h)-10) || ind_sp == (GO(ind_s,ind_phi,ind_h)-30))
                                                    simul(id,t).stat = 1; 
                                                    simul(id,t).eta = opt_eta_O(ind_s,ind_phi,ind_h);
                                                    eta = simul(id,t).eta;
        
                                                    simul(id,t).s_ind = ind_sp;
                                                    ind_s = simul(id,t).s_ind;
                                                    ind_hp = opt_ind_hp_N(ind_s,ind_phi,ind_h);
                                                    hp = h_grid(ind_hp);
                                                    simul(id,t).wage = eta*(VN(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_s,ind_phi,ind_hp)*opt_x_N(ind_s,ind_phi,ind_hp)...
                                                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_s,ind_phi,ind_hp))*(eta*VN(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                                end
                                            end
                                        else                                                          % O - Search OTJ for O
                                            for ind_sp = 1:ns
                                                if (ind_sp== (GO(ind_s,ind_phi,ind_h)-20))
                                                    simul(id,t).stat = 2; 
                                                    simul(id,t).eta = opt_eta_O(ind_s,ind_phi,ind_h);
                                                    eta = simul(id,t).eta;
        
                                                    simul(id,t).s_ind = ind_sp;
                                                    ind_s = simul(id,t).s_ind;
                                                    ind_hp = opt_ind_hp_O(ind_s,ind_phi,ind_h);
                                                    hp = h_grid(ind_hp);
                                                    simul(id,t).wage = eta*(VO(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_s,ind_phi,ind_hp)*opt_x_O(ind_s,ind_phi,ind_hp)...
                                                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_s,ind_phi,ind_hp))*(eta*VO(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp)); 
                                                end
                                            end
                                        end
                                    else                                         % O last period - no separation shock - does NOT successfully searchs OTJ
                                        simul(id,t).stat = simul(id,t-1).stat;
                                        simul(id,t).eta = simul(id,t-1).eta;
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_s,ind_phi,ind_h);  
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_s,ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_s,ind_phi,ind_hp)*opt_x_O(ind_s,ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_s,ind_phi,ind_hp))*(eta*VO(ind_s,ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                                    end
                                end % End if receives separation shock when OL last period
                            end % End if for stat (when not dead or died before period start) last period
                        end         % End if dies at end of last period
                    end    % End if dead in last period
                end % End if t==1
            end % End loop over time periods in panel
        end % End loop over ids in simul_size
        
        
        
        count_O =0.0;
        count_N =0.0;
        
        count_E =0.0;
        count_U =0.0;
        for id=1:simul_size
            for t = 1:nt
                if (simul(id,t).stat==2)
                    count_O = count_O + 1.0;
                elseif (simul(id,t).stat==1)
                    count_N = count_N + 1.0;
                end
                
                if (simul(id,t).stat==0) 
                    count_U = count_U + 1.0;
                elseif (simul(id,t).stat==1 || simul(id,t).stat==2)
                    count_E = count_E + 1.0;
                end
            end
        end
        
        if policy_ind==1
            %% =============================================================================================================== %%
            %========================================== Compute Moments of Interest ============================================%
            %================================================================================================================= %%
                
            % ----------- m1: JTJ Probability
            jtj_rate = 0.0;
            m_sum = 0.0;
            for ind_phi=1:nphi
                for ind_h = 1:nh
                    for ind_s = 1:ns
            
                        ind_hp_N = opt_ind_hp_N(ind_s,ind_phi,ind_h); 
                        ind_hp_O = opt_ind_hp_O(ind_s,ind_phi,ind_h);
                        
                        jtj_rate = jtj_rate + lambda*opt_p_N(ind_s,ind_phi,ind_hp_N)*(1.0-delta)*dist_N(ind_s,ind_phi,ind_h)*(1.0-nu);
                        jtj_rate = jtj_rate + lambda*opt_p_O(ind_s,ind_phi,ind_hp_O)*(1.0-delta)*dist_O(ind_s,ind_phi,ind_h)*(1.0-nu);
            
                        jtj_rate = jtj_rate + opt_p_U(ind_phi,ind_hp_N)*delta*dist_N(ind_s,ind_phi,ind_h)*(1.0-nu);
                        jtj_rate = jtj_rate + opt_p_U(ind_phi,ind_hp_O)*delta*dist_O(ind_s,ind_phi,ind_h)*(1.0-nu);
                    
                        m_sum = m_sum + dist_N(ind_s,ind_phi,ind_h)*(1.0-nu) + dist_O(ind_s,ind_phi,ind_h)*(1.0-nu);
                    end
                end 
            end 
               
            if (m_sum>0.0) 
                jtj_rate = (jtj_rate/m_sum);
            else
                jtj_rate = 999.0;
            end
            m1 = jtj_rate;
            
            % ----------- m2: Unemployment Rate
            m2 = (sum(sum(dist_U(:,:)))/( sum(sum(sum(dist_N(:,:,:)))) + sum(sum(sum(dist_O(:,:,:)))) + sum(sum(dist_U(:,:))) ))*100.0;
            
            % ----------- m3: Total Percent Outsourced
            m3 = (( sum(sum(sum(dist_O(:,:,:)))) )/(sum(sum(sum(dist_N(:,:,:)))) + sum(sum(sum(dist_O(:,:,:))))))*100.0;
            
            % ----------- m4: Ratio Outsourced With to Without Bachelor's Degree
            DO_phi1 = (( sum(sum(dist_O(:,1,:))) )/( sum(sum(dist_O(:,1,:))) + sum(sum(dist_N(:,1,:))) ))*100.0;
            DO_phi2 = (( sum(sum(dist_O(:,2,:))) )/( sum(sum(dist_O(:,2,:))) + sum(sum(dist_N(:,2,:))) ))*100.0;
            m4 = DO_phi1/DO_phi2;
            
            % ----------- m5: Avg. Quarterly Separation Rate
            delta_rate = 0.0;
            m_sum = 0.0;
            
            for ind_phi=1:nphi
                for ind_h = 1:nh
                    for ind_s = 1:ns
                    
                        ind_hp_N = opt_ind_hp_N(ind_s,ind_phi,ind_h); 
                        ind_hp_O = opt_ind_hp_O(ind_s,ind_phi,ind_h);    
            
                        delta_rate = delta_rate + (1.0-opt_p_U(ind_phi,ind_hp_N))*delta*dist_N(ind_s,ind_phi,ind_h)*(1.0-nu);
                        delta_rate = delta_rate + (1.0-opt_p_U(ind_phi,ind_hp_O))*delta*dist_O(ind_s,ind_phi,ind_h)*(1.0-nu);
            
                        m_sum = m_sum + dist_N(ind_s,ind_phi,ind_h)*(1.0-nu) + dist_O(ind_s,ind_phi,ind_h)*(1.0-nu);
                    end
                end 
            end 
            m5 = delta_rate/m_sum;
            
            % ----------- m6: Average Wage Loss After Unemployment Spell   
            wg_uspell = 0;
            m_sum = 0;
            for id = 1:simul_size
                stat_u = 0;
                for t = 1:nt-2 % Need to observe employed wage, unemployment, and then employed wage again
                    
                    if (simul(id,t).stat~=999)
                        if (stat_u==1 && simul(id,t).stat~=0) % Was Unemployed but found Employment
                            w_after = simul(id,t).wage;
                            wg_uspell = wg_uspell + (log(w_after) - log(w_before))*100;
                            m_sum = m_sum +1;
                            stat_u = 0;
                        end
            
                        if (simul(id,t).stat~=0 && simul(id,t+1).stat==0) % Employed in t and unemployed in t+1
                            w_before = simul(id,t).wage;
                            stat_u = 1;
                        end 
                    end % End if has not exited the model (if stat~=999)
                end
            end
            wg_uspell = wg_uspell/m_sum;
            m6 = wg_uspell;
            
            % ----------- m7: Wage Dispersion p90/p50
            count_wages = 0;
            for id = 1:simul_size
                for t = 1:nt       
                    if (simul(id,t).stat~=0 && simul(id,t).stat~=999)
                        count_wages = count_wages + 1;
                    end    
                end
            end
                    
            % --- First save sample of wages
            wage_sim = zeros(1,count_wages);
            count_wages = 0;
            for id = 1:simul_size
                for t = 1:nt       
                    if (simul(id,t).stat~=0 && simul(id,t).stat~=999)
                        count_wages = count_wages + 1;
                        wage_sim(count_wages) = simul(id,t).wage;
                    end    
                end
            end
            p90_wage = prctile(wage_sim, 90);
            p50_wage = prctile(wage_sim, 50);
                    
             m7 = p90_wage/p50_wage;
                    
            
            
            
            
            %----------- m8: Average Annual Real Wage Growth
            % average wage growth over 1 year
            wg_sum2 = 0.0;
            m_sum = 0.0;
            for id = 1:simul_size
                for t = 1:(nt-4)
                    if (simul(id,t).stat~=0 && simul(id,t+4).stat~=0 && simul(id,t).stat~=999 && simul(id,t+4).stat~=999) 
                        w1 = simul(id,t).wage;
                        w5 = simul(id,t+4).wage; 
                        wg_sum2 = wg_sum2 + (log(w5) - log(w1))*100.0;
                        m_sum = m_sum + 1.0;
                    end 
                end 
            end 
            wg_sum2 = wg_sum2/m_sum;
            m8 = wg_sum2;
            
            
            % ----------- m9: College Wage Premium
            avg_wage_phi1 = 0.0;
            avg_wage_phi2 = 0.0;
            sum_phi1= 0.0;
            sum_phi2= 0.0;
            
            for id = 1:simul_size
                for t = 1:nt
                    
                    if (simul(id,t).stat~=0 && simul(id,t).stat~=999)   % If Not Unemployed or dead
                        if (simul(id,t).phi_ind == 1) 
                            avg_wage_phi1 = avg_wage_phi1 + simul(id,t).wage;
                            sum_phi1 = sum_phi1 + 1.0;
                        else
                            avg_wage_phi2 = avg_wage_phi2 + simul(id,t).wage;
                            sum_phi2 = sum_phi2 + 1.0;
                        end 
                    end 
                end 
            end 
            
            avg_wage_phi1 = avg_wage_phi1/sum_phi1;
            avg_wage_phi2 = avg_wage_phi2/sum_phi2;
            
            m9 = avg_wage_phi2/avg_wage_phi1;
            
            
            %---------- m10 - Unemployment Rate Among Age 20
            count_E =0.0;
            count_U =0.0;
            for id = 1:simul_size
                for t = 1:1 % age = 20
                    if (simul(id,t).stat~=999) 
                        
                        if (simul(id,t).stat==0) 
                            count_U = count_U + 1.0;
                        elseif (simul(id,t).stat==1 || simul(id,t).stat==2) 
                            count_E = count_E + 1.0;
                        end
                        
                    end
                end 
            end 
            m10 = (count_U/(count_E + count_U))*100.0;
            
            
            
            %--------------------- m11: Vacancy rate: vacancies/(vacancies + filled jobs)
            % in each submarket - theta = v/u , so v = theta*u
            vacancies = 0.0;
            for ind_phi = 1:nphi   
                for ind_h = 1:nh
            
                 
                    vacancies = vacancies + opt_theta_U(ind_phi,ind_h)*dist_U(ind_phi,ind_h); 
                    for ind_s = 1:ns
                        if (~isreal(opt_theta_N(ind_s,ind_phi,ind_h))) 
                        else
                            vacancies = vacancies + opt_theta_N(ind_s,ind_phi,ind_h)*lambda*dist_N(ind_s,ind_phi,ind_h);
                        end 
                        
                        if (~isreal(opt_theta_O(ind_s,ind_phi,ind_h))) 
                        else
                            vacancies = vacancies + opt_theta_O(ind_s,ind_phi,ind_h)*lambda*dist_O(ind_s,ind_phi,ind_h);
                        end 
                    end
                end 
            end  % end loop over ind_phi
            vacancy_rate = (vacancies/(vacancies + sum(sum(sum(dist_N(:,:,:)))) + sum(sum(sum(dist_O(:,:,:)))) ))*100.0;
            m11 = vacancy_rate;
            
            %----------- m12: New Entrant vs. Older Wage Ratios
            % Look at Wage Ratio 25-34/35-44 to <25/25-34
            avg_lt25_wage = 0;
            count_lt25 =0;
            avg_25_34_wage = 0;
            count_25_34 = 0;
            avg_35_44_wage = 0;
            count_35_44 = 0;
            
            for id = 1:simul_size
                for t = 1:nt
                    if (simul(id,t).stat~=999)
                        if (t < 21) % Age < 25 yrs
                            if (simul(id,t).stat~=0)
                                avg_lt25_wage = avg_lt25_wage + simul(id,t).wage;
                                count_lt25 = count_lt25 + 1;
                            end
            
                        elseif (t>=21 && t <61) % Age 25-34 yrs
                            if (simul(id,t).stat~=0)
                                avg_25_34_wage = avg_25_34_wage + simul(id,t).wage;
                                count_25_34 = count_25_34 + 1;
                            end
            
                        elseif (t>=61 && t <101) % Age 35-44 yrs
                            if (simul(id,t).stat~=0)
                                avg_35_44_wage = avg_35_44_wage + simul(id,t).wage;
                                count_35_44 = count_35_44 + 1;
                            end
                        else
                        end
                    end
                end
            end
            avg_lt25_wage = avg_lt25_wage/count_lt25;
            avg_25_34_wage = avg_25_34_wage/count_25_34;
            avg_35_44_wage = avg_35_44_wage/count_35_44;
            m12 = (avg_35_44_wage/avg_25_34_wage)/(avg_25_34_wage/avg_lt25_wage);
            
            mm_D4 = [m1;m2;m3;m4;m11;m5;m6;m7;m8;m12;m9;m10];
            param_D4 = [lambda;b;muA;c(1);c(2);delta;d;chi;alpha(1);alpha(2);z(2);h_sigma];
                 
            %-------------- Model Selection into Outsourcing --------------------%
            perc_college_N = (( sum(sum(dist_N(:,2,:))))/( sum(sum(sum(dist_N(:,:,:)))) ))*100;
            perc_college_O = (( sum(sum(dist_O(:,2,:))))/( sum(sum(sum(dist_O(:,:,:)))) ))*100;
                    
            % Selection: Average h
            avg_h_N = 0;
            avg_h_O = 0;
            for ind_h = 1:nh
                h = h_grid(ind_h);
                    
                avg_h_N = avg_h_N + h*(sum(sum(dist_N(:,:,ind_h))) );
                avg_h_O = avg_h_O + h*(sum(sum(dist_O(:,:,ind_h))) );
            end
            avg_h_N = avg_h_N/( sum(sum(sum(dist_N(:,:,:)))) );
            avg_h_O = avg_h_O/( sum(sum(sum(dist_O(:,:,:)))) );
                    
            % Alpha
            alpha_O = 0.0;
            m_sum_O = 0.0;
            alpha_N = 0.0;
            m_sum_N = 0.0;
            for ind_phi=1:nphi
                for ind_h = 1:nh
                    for ind_s = 1:ns
                    
                        alpha_N = alpha_N + alpha(ind_s)*dist_N(ind_s,ind_phi,ind_h)*(1.0-nu);
                        m_sum_N = m_sum_N + dist_N(ind_s,ind_phi,ind_h)*(1.0-nu); 
                                
                        alpha_O = alpha_O + alpha(ind_s)*dist_O(ind_s,ind_phi,ind_h)*(1.0-nu);
                        m_sum_O = m_sum_O + dist_O(ind_s,ind_phi,ind_h)*(1.0-nu);
                    end
                end 
            end 
            alpha_O = (alpha_O/m_sum_O);
            alpha_N = (alpha_N/m_sum_N);
        
            %------------------- Optimal Investment Rule: h'-h --------------------%
            %--- Want to plot Investment Difference: Non-Outsourced, Non-college, High versus low separation types %
            %--------------------------------------: Non-college, low separation, Outsourced vs. Non-outsourced    %
                    
            opt_hp_N = h_grid(opt_ind_hp_N);
            opt_hp_O = h_grid(opt_ind_hp_O);
                    
            h_diff_L(:) = reshape(opt_hp_N(1,1,:),[1,nh]) - h_grid(:)';
            h_diff_H(:) = reshape(opt_hp_N(2,1,:),[1,nh]) - h_grid(:)';
                    
            h_diff_N(:) = reshape(opt_hp_N(1,1,:),[1,nh]) - h_grid(:)';
            h_diff_O(:) = reshape(opt_hp_O(1,1,:),[1,nh]) - h_grid(:)';
                    
            avg_diff_L(1) = mean(h_diff_L(1:108));        %0-0.25
            avg_diff_L(2) = mean(h_diff_L(109:215));      %0.25-0.5
            avg_diff_L(3) = mean(h_diff_L(216:322));      %0.5-0.75
            avg_diff_L(4) = mean(h_diff_L(323:429));      %0.75-1
                        
            avg_diff_L(5) = mean(h_diff_L(430:536));      %1-1.25
            avg_diff_L(6) = mean(h_diff_L(537:643));      %1.25-1.5
            avg_diff_L(7) = mean(h_diff_L(644:750));      %1.5-1.75
            avg_diff_L(8) = mean(h_diff_L(751:857));      %1.75-2
                        
            avg_diff_L(9) = mean(h_diff_L(858:964));      %2-2.25
            avg_diff_L(10) = mean(h_diff_L(965:1072));    %2.25-2.5
            avg_diff_L(11) = mean(h_diff_L(1073:1179));   %2.5-2.75
            avg_diff_L(12) = mean(h_diff_L(1180:1286));   %2.75-3
                        
            avg_diff_H(1) = mean(h_diff_H(1:108));        %0-0.25
            avg_diff_H(2) = mean(h_diff_H(109:215));      %0.25-0.5
            avg_diff_H(3) = mean(h_diff_H(216:322));      %0.5-0.75
            avg_diff_H(4) = mean(h_diff_H(323:429));      %0.75-1
                        
            avg_diff_H(5) = mean(h_diff_H(430:536));      %1-1.25
            avg_diff_H(6) = mean(h_diff_H(537:643));      %1.25-1.5
            avg_diff_H(7) = mean(h_diff_H(644:750));      %1.5-1.75
            avg_diff_H(8) = mean(h_diff_H(751:857));      %1.75-2
                        
            avg_diff_H(9) = mean(h_diff_H(858:964));      %2-2.25
            avg_diff_H(10) = mean(h_diff_H(965:1072));    %2.25-2.5
            avg_diff_H(11) = mean(h_diff_H(1073:1179));   %2.5-2.75
            avg_diff_H(12) = mean(h_diff_H(1180:1286));   %2.75-3
                      
            avg_diff_N(1) = mean(h_diff_N(1:108));        %0-0.25
            avg_diff_N(2) = mean(h_diff_N(109:215));      %0.25-0.5
            avg_diff_N(3) = mean(h_diff_N(216:322));      %0.5-0.75
            avg_diff_N(4) = mean(h_diff_N(323:429));      %0.75-1
                        
            avg_diff_N(5) = mean(h_diff_N(430:536));      %1-1.25
            avg_diff_N(6) = mean(h_diff_N(537:643));      %1.25-1.5
            avg_diff_N(7) = mean(h_diff_N(644:750));      %1.5-1.75
            avg_diff_N(8) = mean(h_diff_N(751:857));      %1.75-2
                        
            avg_diff_N(9) = mean(h_diff_N(858:964));      %2-2.25
            avg_diff_N(10) = mean(h_diff_N(965:1072));    %2.25-2.5
            avg_diff_N(11) = mean(h_diff_N(1073:1179));   %2.5-2.75
            avg_diff_N(12) = mean(h_diff_N(1180:1286));   %2.75-3
                        
            avg_diff_O(1) = mean(h_diff_O(1:108));        %0-0.25
            avg_diff_O(2) = mean(h_diff_O(109:215));      %0.25-0.5
            avg_diff_O(3) = mean(h_diff_O(216:322));      %0.5-0.75
            avg_diff_O(4) = mean(h_diff_O(323:429));      %0.75-1
                        
            avg_diff_O(5) = mean(h_diff_O(430:536));      %1-1.25
            avg_diff_O(6) = mean(h_diff_O(537:643));      %1.25-1.5
            avg_diff_O(7) = mean(h_diff_O(644:750));      %1.5-1.75
            avg_diff_O(8) = mean(h_diff_O(751:857));      %1.75-2
                        
            avg_diff_O(9) = mean(h_diff_O(858:964));      %2-2.25
            avg_diff_O(10) = mean(h_diff_O(965:1072));    %2.25-2.5
            avg_diff_O(11) = mean(h_diff_O(1073:1179));   %2.5-2.75
            avg_diff_O(12) = mean(h_diff_O(1180:1286));   %2.75-3

        end
        %=============== Estimate Aggregate Effects of Removing Outsourcing Option =================%
        % Compute and save moments of interest for comparison
        %-- Unemployment rate
        policy_urate(policy_ind) = (sum(sum(dist_U(:,:)))/( sum(sum(sum(dist_N(:,:,:)))) + sum(sum(sum(dist_O(:,:,:)))) + sum(sum(dist_U(:,:))) ))*100;
        
        %----------- Output from Employment (GDP)
        policy_GDP(policy_ind) = 0;
        for ind_s = 1:ns
            for ind_phi = 1:nphi
                for ind_h = 1:nh
                    h = h_grid(ind_h);
                    policy_GDP(policy_ind) = policy_GDP(policy_ind) + (z(ind_phi)*h^alpha(ind_s))*(dist_N(ind_s,ind_phi,ind_h) + dist_O(ind_s,ind_phi,ind_h));
                end
            end
        end
        
        %---------- Average Human Capital
        h_sum_total = 0;
        policy_avg_hc(policy_ind) = 0;
        for ind_s = 1:ns
            for ind_phi = 1:nphi
                for ind_h = 1:nh
                    h = h_grid(ind_h);
        
                    policy_avg_hc(policy_ind) = policy_avg_hc(policy_ind) + h*(dist_U(ind_phi,ind_h) + dist_N(ind_s,ind_phi,ind_h) + dist_O(ind_s,ind_phi,ind_h));
                    h_sum_total = h_sum_total + dist_U(ind_phi,ind_h) + dist_N(ind_s,ind_phi,ind_h) + dist_O(ind_s,ind_phi,ind_h);
                end
            end
        end
        policy_avg_hc(policy_ind) = policy_avg_hc(policy_ind)/h_sum_total;
        
        %-------- Average Wage
        count_wage_total = 0;
        policy_avg_wage(policy_ind) = 0;
        for id = 1:simul_size
            for t = 1:nt
                if (simul(id,t).stat~=999 && simul(id,t).stat~=0) % If employed
                    policy_avg_wage(policy_ind) = policy_avg_wage(policy_ind) + simul(id,t).wage;
                    count_wage_total = count_wage_total + 1;
        
                end
            end
        end
        policy_avg_wage(policy_ind) = policy_avg_wage(policy_ind)/count_wage_total;
            
        %--------- Total Utility
        util = 0;
        for ind_s = 1:ns
            for ind_phi = 1:nphi
                for ind_h = 1:nh
                    h = h_grid(ind_h);
                    hp_N = h_grid(opt_ind_hp_N(ind_s,ind_phi,ind_h));
                    hp_O = h_grid(opt_ind_hp_O(ind_s,ind_phi,ind_h));
                    % Compute production minus investment cost
                    util = util + dist_U(ind_phi,ind_h)*(b*z(ind_phi)*(h^alpha(ind_s)))...
                        + z(ind_phi)*(h^alpha(ind_s))*(dist_N(ind_s,ind_phi,ind_h) + dist_O(ind_s,ind_phi,ind_h))...
                        - (((hp_N - h*(1-d))*chi)^gamma)*dist_N(ind_s,ind_phi,ind_h) - (((hp_O - h*(1-d))*chi)^gamma)*dist_O(ind_s,ind_phi,ind_h);
                    % Now subtract out vacancy posting costs
                    % Now subtract out vacancy posting costs
                    if (GU(ind_phi,ind_h)<20 || GU(ind_phi,ind_h)>=30) % Search in unemployment for non-outsourced job
                        if (GU(ind_phi,ind_h)<20)
                            ind_s_search = GU(ind_phi,ind_h)-10;
                        else
                            ind_s_search = GU(ind_phi,ind_h)-30;
                        end
                        if (~isnan(opt_theta_U(ind_phi,ind_h)) && isreal(opt_theta_U(ind_phi,ind_h)))
                            util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*c(ind_s_search);
                        end
                    else                                              % Search in unemployment for outsourced job
                        ind_s_search = GU(ind_phi,ind_h)-20;
                        if (~isnan(opt_theta_U(ind_phi,ind_h)) && isreal(opt_theta_U(ind_phi,ind_h)))
                            util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*(c(ind_s_search) + A_staffing(ind_phi,ind_h));
                        end
                    end
    
                    if (GN(ind_s,ind_phi,ind_h)<20 || GN(ind_s,ind_phi,ind_h)>=30) % Search from non-outsourced job for non-outsourced job
                        if (GN(ind_s,ind_phi,ind_h)<20)
                            ind_s_search = GN(ind_s,ind_phi,ind_h)-10;
                        else
                            ind_s_search = GN(ind_s,ind_phi,ind_h)-30;
                        end
    
                        if (~isnan(opt_theta_N(ind_s,ind_phi,ind_h)) && isreal(opt_theta_N(ind_s,ind_phi,ind_h)))
                            util = util - dist_N(ind_s,ind_phi,ind_h)*opt_theta_N(ind_s,ind_phi,ind_h)*lambda*c(ind_s_search);
                        end
                    else                                              % Search from outsourced job for non-outsourced job
                        ind_s_search = GN(ind_s,ind_phi,ind_h)-20;
                        if (~isnan(opt_theta_N(ind_s,ind_phi,ind_h))&& isreal(opt_theta_N(ind_s,ind_phi,ind_h)))
                            util = util - dist_N(ind_s,ind_phi,ind_h)*opt_theta_N(ind_s,ind_phi,ind_h)*lambda*(c(ind_s_search) + A_staffing(ind_s,ind_phi,ind_h));
                        end
                    end
    
                    if (GO(ind_s,ind_phi,ind_h)<20 || GO(ind_s,ind_phi,ind_h)>=30) % Search from outsourced job for non-outsourced job
                        if (GO(ind_s,ind_phi,ind_h)<20)
                            ind_s_search = GO(ind_s,ind_phi,ind_h)-10;
                        else
                            ind_s_search = GO(ind_s,ind_phi,ind_h)-30;
                        end
                        if (~isnan(opt_theta_O(ind_s,ind_phi,ind_h)) && isreal(opt_theta_O(ind_s,ind_phi,ind_h)))
                            util = util - dist_O(ind_s,ind_phi,ind_h)*opt_theta_O(ind_s,ind_phi,ind_h)*lambda*c(ind_s_search);
                        end
                    else
                        ind_s_search = GO(ind_s,ind_phi,ind_h)-20;
                        if (~isnan(opt_theta_O(ind_s,ind_phi,ind_h)) && isreal(opt_theta_O(ind_s,ind_phi,ind_h)))
                            util = util - dist_O(ind_s,ind_phi,ind_h)*opt_theta_O(ind_s,ind_phi,ind_h)*lambda*(c(ind_s_search) + A_staffing(ind_s,ind_phi,ind_h));
                        end
                    end
                end
            end
        end
        policy_util(policy_ind) = util;
        
        %---------- New Entrant Value
        policy_entrant(policy_ind) = pr_phi(1)*U_hat(1,1) + pr_phi(2)*U_hat(2,1);
    
    end

    ppc_urate = policy_urate(2)-policy_urate(1);
    pc_GDP = ((policy_GDP(2) - policy_GDP(1))/policy_GDP(1))*100;
    pc_avg_hc = ((policy_avg_hc(2) - policy_avg_hc(1))/policy_avg_hc(1))*100;
    pc_avg_wage = ((policy_avg_wage(2) - policy_avg_wage(1))/policy_avg_wage(1))*100;
    pc_total_util = ((policy_util(2) - policy_util(1))/policy_util(1))*100;
    pc_entrant = ((policy_entrant(2) - policy_entrant(1))/policy_entrant(1))*100;
    
    policy_D4 = [ppc_urate;pc_GDP;pc_avg_hc;pc_avg_wage;pc_total_util;pc_entrant];

end % end function
